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Abstract. The Robinson-Foulds (RF) distance is by far the most widely used measure of 
dissimilarity between trees. Although the distribution of these distances has been investigated 
for twenty years, an algorithm that is explicitly polynomial time has yet to be described for 
computing this distribution (which is also the distribution of trees around a given tree under 
the popular Robinson-Foulds metric). In this paper we derive a polynomial-time algorithm for 
this distribution. We show how the distribution can be approximated by a Poisson distribution 
determined by the proportion of leaves that lie in 'cherries' of the given tree. We also describe 
how our results can be used to derive normalization constants that are required in a recently- 
proposed maximum likelihood approach to supertree construction. 



1. Introduction 

Tree comparison metrics are widely used in phylogenetics for comparing evolutionary trees [21 [5] 
and for performing statitistical tests - for example, to test whether two trees are 'significantly 
different' from each other than one might expect if one or both trees were randomly chosen OH]. 
In order to address these statistical questions one needs to determine the distribution of the 
metric under some null model (see, for example, [5lI6])- The symmetric difference or Robinson- 
Foulds metric is the most widely used measure of differences between phylogenetic trees, and its 
distribution is particularly attractive to study. In a landmark paper [S], the authors described 
this distribution of trees relative to a fixed reference tree via a system of generating functions. 
This allowed the authors to calculate the distribution explicitly for small trees and provided a 
tool for analytic results on this distribution in later work by others. 

However, the approach described in p51 does not immediately appear to provide a polynomial- 
time algorithm for computing this distribution, and for larger trees their approach may be com- 
putationally prohibitive. In this paper, we describe how to calculate the distribution of the 
Robinson-Foulds metric relative to a fixed tree. We also show how the distribution can be ap- 
proximated by a Poisson distribution whose parameter depends on just one aspect of tree shape 
- the number of 'cherries'. 
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Our investigation into the distribution of the metric has also been motivated by its relevance 
to a recent approach for 'supertree' construction that is based on maximum likelihood [S]- In 
particular, our algorithm allows the normalization constants in the likelihood calculations to be 
computed explicitly. We describe how these normalization constants depend weakly on aspects 
of the shape of the tree - for example, how many 'cherries' the tree has. We start by recalling 
some terminology. 

Let X be a finite set. A phylogenetic tree with leaf set X is a tree with its degree one vertices 
(leaves) labelled bijectively by elements of X and whose remaining vertices have degree at least 
three. We use V{T) and E{T) to denote the set of nodes (vertices) and edges of T. Let V(T) 
denote the set of internal (non-leaf) nodes of T and let E{T) be the set of edges in E{T) that 
have both endpoints in V{T), the internal edges. 

A phylogenetic tree is fully resolved if every internal vertex has degree three. Following [3] we let 
PT{n) denote the set of phylogenetic trees on the finite set X = {1, 2, . . . , n} and BPT{n) the 
set of fully resolved ('binary') trees in PT{n) (two trees in BPT{6) are shown in Fig. 1). The 
number of trees in BPT{n) is denoted b{n) and is given by: 

n 

(1) 6(n) (2n-5)!! = J|(2fc-5) n > 3, 

fc=3 

see [H]. For convenience, we let /3(m) denote the number of fully resolved trees with exactly m 
internal edges, so: 

m+3 

(2) /3(m) ^ 6(m + 3) = J]^ (2fc - 5) m > 0. 

fc=3 

Every edge e G E{T) induces a bipartition or split of the leaf set X corresponding to the labels 
present in the two connected components remaining when the edge e is removed. Let 7r(T, e) 
denote this bipartition, which we consider unordered. We let c(T) denote the set of all bipartitions 
obtained by removing different edges of e. Hence |c(r)| < 2n — 3, the maximum number of edges 
in a phylogenetic tree, and \c{T)\ — 2n — S exactly when T is fully resolved. A bipartition is 
trivial if it separates a single element from all other elements; trivial bipartitions correspond to 
the edges in the tree that are external, meaning that they are incident with a leaf of the tree. 
A cherry of a fully resolved phylogenetic tree T is a pair of leaves that forms one half of a split 
of T (i.e. a pair of leaves whose incident edges contain a common vertex). In Fig. 1 the pairs 
(1,2) and (5,6) form cherries in both trees, while the right-hand tree has an additional cherry 
(3,4). 

The symmetric difference metric is defined on PT(n), and hence on BPT(n), by: 

d{TuT2) = |c(Ti) Ac(r2)|. 
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Note that this number is always even when Ti and T2 are both in BPT{n), since, for any two 
trees in PT{n), we have d{Ti,T2) = \c{Ti)\ + \c{T2)\ - 2|c(Ti) n c(T2)|, and if Ti,T2 G BPT{n) 
then |c(ri)| — 10(12)1 = 2n — 3. As an example, the two trees shown in Fig. 1 have a distance 
value of 2 since the splits {1, 2, 3}|{4, 5, 6} and {3, 4}|{1, 2, 5, 6} each occur in just one tree. 

The metric was introduced by Bourque [T] and generalised by Robinson and Foulds [7]. As all 
phylogenetic trees contain all trivial splits, the maximum possible distance between two trees is 
2(n — 3), which is twice the maximum number of internal edges. 



2. Computing the distribution of the Robinson-Foulds metric 

For each T G PT{n), let 6,„(r) denote the number of trees T' G BPTin) for which d(T, T') = m. 
As d is a metric, 60 (r) = 1. A recursive formula for the generating function of hm{T) is given 
in [3] and [llj . As far as we could deduce, the formula does not provide a polynomial time 
algorithm for computing the bm{T) values, due to an exponential explosion in the number of 
subcases. 

Instead we use an alternative approach, applying results of llj. Let qs{T) denote the number of 
trees in BPT{n) that share exactly s internal splits with T. Then for all m = 0, 2, 4, ... , 2(n — 3), 
we have: 

(3) hra{T) = qn-3-7n/2{T). 

Define the polynomial 

n-3 



(4) q{T,x) = Y,<lAT)a 



s=0 

Let E C E{T) denote a subset of the set of internal edges of T. The forest T — E has exactly 
\E\ + 1 components Fi, F2, . . . , F\e\+i- We use E{Fi) as a short-hand for the edges of E{T) that 
are contained in Fi. 

Define 

(5) ne{t)= n mm)\) 

1=1 

Then Ne{T) equals the quantity defined in [TT] (here assuming that T is fully resolved). 

For s > define 

ECE(T) 
\E\=s 
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the sum of Ne over all subsets E C E{T) of cardinality s. For example, 7'o(T) equals f3{\E{T)\) = 
(3{n — 3). It was shown in [TT] that the generating function 

i?(T,x)=^r,(r)a;^ 

satisfies the identity 

(6) q{T,x)^R{T,x~l). 

In what follows we derive a formula to evaluate the coefficients rs{T) so that we can compute 
the coefficients b„i{T) via ^ and 

As usual, the computation applies dynamic programming. Let vq be the node adjacent to leaf 
n. Delete leaf n and make vq the root of the tree, so that now every internal node has exactly 
two children. For each internal node v let T„ denote the subtree of T containing v and all of its 
descendants. Given a subset E C E{Ty), we define NsiTy) as in where Fi, . . . , F\e\+i will 
now be components oi Ty — E instead oi T — E. We let k{v, E) denote the number of edges in 
the component oiT^ — E containing v. For s, fc > 0, we let £{v, s, k) denote the set of all subsets 
E C E{Ty) such that \E\ = s and k{v, E) = k. Define 

(7) R{v,s,k)^ J2 ^Ein) 

Ee£{v,s,k) 

so that if Vq is the root of T and s > 0, we have: 

(8) r,(T) = ^i?(z;o,s,fc). 

We now derive a recursion for R{v, s, fc). As is customary, an empty summation equals zero. 



Lemma 1. Suppose that v £ V{T). Then 
(9) R{v,0,k) = 



m tfk=\E{n)\; 
otherwise. 



Lemma 2. Suppose that s > 1. For all v £ E{T) let riy = \E{Ty)\. 

(1) If k > Uy then R{v, s, fc) = 0. 

(2) If V £ V{T) has no children in V and s > 1 then R{v, s, fc) = 0. 

(3) If V £ V{T) has one child vi in V then 

s, fc — l)(2fc + 1) otherwise; 



(10) R{v,s,k) 
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(4) IfvG V{T) has two children vi,V2 in V{T) then 



s-2 



(11) R{v, s,0)=^ {J2 Ri->Ji^si,ki) ^(^2, s - 2 - Si, fca) 

si=o yfei>o / \fe2>o 

(5) Ifv€ V{T) has two children vi,V2 in V{T) and k >1 then 

R{v,s,k) = ^ I ^ i?(i;i,si,A:i) 1 i?(^;2,s-l-si,fc-l)/3(A;)//3(fc-l) 

si=o yfci>o / 

+ fc f E R{v2,S2,k2)] R{Vi,S-l-S2,k-l)(}{k)/P{k-l) 

S2=a \k2>0 J 

+ J2 51 -^(^i'^i'^i)-^(^2,s-si,fc-2-A;i^ 



p{k,mk-2-k. 



(12) 



Proof. Parts (1) and (2) follow from the definition of R. 

(3) Let e be the edge from vi to v. When A: = it holds that E G £{v,s,k) if and only if 
E = EiL) {e} for some Ei e £{vi,s — 1, fci), where ki ranges from to s — 1. This gives 
the first case. When A; > 1, the edge e connecting v and vi is absent from every set in 
£{v, s, k). Thus E e £{v, s, k) if and only if £ e £{vi,s, k — 1). 

= NE'{n,){2k + i). 

(4) Let ei, 62 be the edges from v to vi,V2 respectively. Since k = 0, for all E e £{v, s, k), we 
have ei ^ E and €2 ^ E. Thus E e £{v, s, k) if and only if there exists E\ e £{v\,s\, ki) 
and E2 e £{v2,s — 51,^2) for some si, ki,k2 > such that E = Eid E2. For each such 
set E, we have: 7Vb(T^) = Ne^{T^^)Ne2{Tv2)- 

(5) Again, let 61,62 be the edges from v to v\,V2 respectively. For each E G £{v,s,k) with 
fc > 0, exactly one of the following cases holds: 

Case 1: 61 & E but 62 E. This case applies if and only there exists Ei G 
£{vi,si,k\) and E2 G ^(1)2,3 — 1 — si,fc — 1) for some si,fci > such that E = 
EiU E2U {ei}. For such a set E we have 

NEin) = NEATv,)NEM^^^y 



Case 2: ei ^ E but 62 G E. Identical to Case 1 with vi and V2 switched. 
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Case 3: ei G E and 62 G E. This case applies if and only there exists Ei G 
£{vi,si,ki) and E2 G £(w2, s - si, A: - fci - 2) such that E — EiU E2U {ei, 62}. For 
each such set E we have: 

m 



(3{h)(3{k-2-h)- 

a 



Theorem 3. Given a fully resolved tree T on n leaves the coefficients hm{T) can he computed 
in 0(n^) time. 

Proof. Consider a vertex v G V{T). If v has one child in V{T) then we evaluate (fTO|) for all 
s, fc < n — 3 in 0(n^) time. If v has two children in V{T) then we evaluate (fT^ in 0{n'^) time. 

Hence computing all the coefficients ^^(r) takes O(ri^) time. From ([6]), we obtain: 

(13) g™(r) = ^ r r,(r)(-ir™, 



s—ra 



from which we compute the values bm{T) — qn-3-m/2{T). □ 



3. POISSON APPROXIMATION 

When n is large we can approximate the QsiT) values by a Poisson distribution with mean 
Xt := cr/'in where ct denotes the number of cherries of T (recall that a cherry is a pair of 
leaves whose incident edges contain a common vertex). More precisely, we have the following 
result. 

Theorem 4. For any tree T G BPT{n), let Yt he a Poisson random variahle with mean At- 
Then the distributions qs{T) /b{n) (the proportion of trees in BPT{n) that share s nontrivial 
splits with T) and Yt have variational distance that converges to zero as n ^ 00. In particular, 

J2 \qs{T)/b{n) - e-^-A^/,s!| = 0{n-\ 

s>0 

Proof. Let Xt denote the random variable which counts the number of non-trivial splits that 
T shares with a tree T' selected uniformly at random from BPT{n). Thus, ^{Xt ^ s) — 
qs{T)/b{n). Let X^ be defined in the same ways as for Xt but counting only splits that divide 
the leaf set into subsets of size 2 and n — 2. Clearly, X^ < Xt. Moreover, the probability that 
T' shares a split with T that is not of the type counted by X!p is bounded above by a term of 
order n^^ and so we have: 

(14) ¥{Xt X^) ^0{n-^). 
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Now, for any two discrete random variables X and X' an elementary probability argument shows 
that Y.S = - = s)\< 2P(^ / X'), and so: 



^ |P(Xt = s) - P(X^ = s)| < 2P(Xt 7^ X't). 



s>0 



(15) 

Combining p4p and (fT5|) gives: 

(16) ^ |P(Xt s) - P(X^ = s)| = O(n-i). 

s>0 

By the triangle inequality, 

(17) ^ |P(Xt = s)-¥{Yt = s)| < ^ |P(Xt = s)-~V{X't = s)|+Xl |1P(^t = s)-¥{Yt = s)| 

s>0 s>0 s>0 

which, combined with p6p . gives: 

(18) = s) - P(rT = s)l < ^ |P(^T = - P(>T = s)| + O(n-i). 

s>0 s>0 

Thus, to establish Theorem [Hit suffices to show that 



(19) Y\^{X'j, = s)-f{YT = s)\=0{n-^) 

s>0 

Now, by Lemma 3 of [TTl, we have: 

(20) nx'T^s) = Y,{-iy 

r—s 

Furthermore, letting A denote At for brevity, we have 



r\ / ct\ h{n — r) 
sj \ r J b{n) 



»(y^^.) = e-^A7s! = f;(-l)'-+^Q 



A"^ 



Substituting this and ([20)1 into the left-hand side of gives the expression: 

^cyA h{n — r) y 



(21) 



E 



r+s 



r J b{n) 



which, after some algebra, and moving the absolute value inside the second summation, is 
bounded above by: 



(22) 
where 



1 °° 1 



s! ^-^ (r — s] 

s>0 r=s ^ ' 



n[=i(i -vct) 



n;=i(l-(2j + 3)/2n) 
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Using the fact that Ct < ^1/2, and a somewhat tedious case analysis, it can be shown that 
f{n,r) < C/n for a constant C that is independent of r,n. It foUows that 

00 

which estabhshes ([T9|) and thereby the theorem. □ 

Remark If T is selected uniformly at random from BPT{n), then At converges in probability to 
I (since the variance of At is 0{n~^) by Theorem 4(b) of [1]). Thus, TheoremlUcan be viewed 
as a refinement of the main result from [llj that for two trees selected uniformly at random from 
BPT{n) the number of non-trivial splits they share is asymptotically Poisson distributed with 
mean -j. 



(24) Zt. = = ^ 



Application to Likelihood based supertrees 

Rodrigo and Steel [9] recently presented a likelihood framework for constructing consensus trees 
and supertrees. Let LiTi) denote the set of leaves of a (fully resolved) gene tree T^. The 
probability of observing Ti with leaf set C{Ti) = Xi given an estimated species tree or supertree 
T has the form 

(23) Pt.x. m) = Pt(T,) = .^—e-^^'^^"^\^^^^ 

^T\C{Ti) 

where T\C{Ti) denotes the restriction of T to the leaf set T^, and where Pi is a positive constant. 
The normalising constant 

^-ftd(T',T|£Ti) 
T':C(T')=C(Ti) 

is required so that the PtCTi) values sum to 1 over all choices of T,;. One complication with 
this approach is that the normalising functions Zt^ depend on T (more precisely, although Zt^ 
does not depend on how the leaves of T are labeled, it may depend on the shape of T), meaning 
that the constant needs to be computed in order to compare the likelihood values of two trees. 
This was overlooked in 9J, in particular Proposition 1 of that paper may only hold in certain 
cases (for example, if the sets Xi are of size at most 5, or if the j3i values are sufficiently large). 
However, Proposition 1 of [9] can be corrected by replacing the term 

k 

Y,|3d{T^,T\Xi) 

in the statement of that Proposition by 

k 

Y,M{T.,T\X,)+j,{T), 
1=1 
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where 

k 

7,(T) = 5]log(^Tj = log(l + J2 e-''""",„(T)), 

i=l m>0 

and where nm{T) is the number of fuUy resolved phylogenetic trees on leaf set Xi that have 
distance ni from T\Xi. 

In general, normalising constants are difficult to evaluate. When d is the Robinson-Foulds dis- 
tance, however, computing the constant is straight-forward. Suppose that \Xi\ = n and that 
bm{T) has been computed for all m. Then (suppressing the index i) we have: 

T'^BPT(n) 

= 5]6™(r)e-^™. 

m 

which can be evaluated directly from the b„i{T) values, and thereby in polynomial time overall 
in n. 

It is instructive to estimate Zt in two limiting cases - firstly for values of j3 that are close to 0, 
and for values of /3 that are large. In both cases we find that the dominant aspect of the shape of 
T affecting Zt is the number ct of cherries that T has. The experimental performance of these 
approximations is evaluated in the final section. 



3.1. Small values of (5. For (3 close to 0, we exploit the identity e /3™ = i _ pffi + O{0^) and 
write: 

(25) ZT = 6(n)-/?^m6„(r)-fO(/32). 

m 

Now, the first term in in bin) times the expected RF distance (denoted v{Ty) from T to a 
tree that is sleeted uniformly at random from BPTin). From [11] (p. 550), we have: 

(26) .{T)/h[n) = 2n - 6 - 2 5: + ^^^^^ " * + 



6(n) 

where rii is the number of interior edges of T for which the smaller subtree in T — e contains i 
leaves of T . For example, consider the 105 fully resolved trees with six leaves, each of which has 
one of two possible shapes, depending on whether it has two or three cherries (as shown in Fig. 
1). For any such tree Ti with two cherries we have: 

KT2)/6(r.) = 6-|, 
while for any tree Ta with three cherries we have: 

30 

Kr3)/6(n) = 6--. 
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Returning to the general setting, we can expand (|26p and write: 

(27) HT)/Kn) - 2n - 6 - 2^ - 6^^^^^^^^ - O(n-), 

where tt is the number of edges e of T for which one of the subtree of T — e has exactly three 
leaves of T. Notice that: 

Tt < CT, 

since any 3-leaf subtree necessarily contains a cherry; therefore a corollary of is 

(28) iy(T)/b(n) ^2n^6-2—^ 0(n-'^), 

2n — 5 

and so, from ([25]), we have 

Zt = b{n) (l - (3{2n - 6 - 2-^ - 0{n-^)) + 0{P^)) . 

Thus, as /? converges to 0, Zt converges to a constant, and when /3 is close to 0, the small 
difference from this constant is dominated by ct- 

3.2. Large values of /3. When fi is large, let e = e^^'^. Then, 

Zt = 1 + h2(T)e + h^(T)^ + 0{t^). 
Now, 62(T) = 2(n - 3), and from Theorem 2.26 of [10^ we have: 

&4(T)=4(^"2^^ +6(7^-6 + ct). 
Thus if we let ^„,e := 1 + (2n - 3)e + 2{r? - An - 6)e^ then 

Zt-A„,, + 6cTe' + 0(e'). 

Once again we see that in the limit (in this P tends to infinity) Zt converges to a constant, 

and for large values of f3, the small difference from this constant is dominated by ct- 

4. Experimental results 

To study general features of the distribution, and examine the accuracy of the above approx- 
imations, we generated random trees and computed the distribution of the Robinson Foulds 
distance for each tree. The trees were drawn from a uniform distribution, with the number of 
taxa varying from 5 to 50. One thousand replicates were performed for each number of taxa. We 
also constructed an unrooted caterpillar tree and a balanced unrooted tree for every set of taxa. 
A balanced unrooted tree is one that minimises the length of the longest path between any two 
leaves, an example being the right-hand tree in Fig. 1. 

As predicted from the Poisson approximation, the distributions of Robinson-Foulds distances 
from a fixed tree were highly peaked. For all of the trees examined, at least 99% of trees are 
either at distance 2{n — 3), the maximum possible, or distance 2(n — 4). 
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For T G BPT{n), let Nk{T) denote the number of trees in BPT{n) within Robinson- Foulds 
distance k of T; that is, 

k 

^k{T) = b,n{T). 

Then N2{T) — 2(n — 3) + 1, the number of trees that share all but one split with T, together 
with the tree T itself. When k > 2, the value of Nk{T) varies with the shape of T. We observed 
that for all k, Nk{T) was minimised when T is a caterpillar. At the other extreme, Nk{T) was 
almost always maximised when T was balanced, the exception being when T was balanced but 
did not have the maximum number of cherries. 

For each tree, and a range of different values for /?, we computed the exact normalising constant 
Zt- Fig. 2 illustrates the variation in Zt over different values of /?, displayed on a log-log plot. 
The central curve gives the average Zt values for 1000 fifty-taxa trees drawn from a uniform 
distribution, as a function of /3. The small-/? and large-/3 approximate values for Zt are also 
plotted. 

As a function of (3, the normalising constant has three distinct phases. For < (3 < 0.03 the 
normalising constant Zt is close to the total number of fully resolved trees, and is fit well by the 
small-/3 approximation. For /3 > 3.0 the normalising constant approaches 1, and is close to the 
large-/? approximation. Between 0.03 and 3, the Zt value changes quickly as a function of /?. In 
this interval, neither of the above approximations work well. 

As we observed above, to correctly compute the likelihood for a supertree under the model of 
[S] we need to compute Zt for every distinct supertree T. Even though this calculation take 
polynomial time, it is still extremely expensive computationally, particularly considering that 
millions of candidate supertrees may be considered. We ask, then, the extent to which this 
computation is strictly necessary. In particular, if we ignore the normalising constant when 
comparing likelihoods, would the relative likelihood ordering of distinct trees change. The key 
question is then to determine how much the normalisation constants Zt vary. If the difference 
is sufficiently small then there will be no impact of ignoring the differences between normalising 
constants. 

For a given value of (3 define the range of Zt to be the ratio of the largest to the smallest Zt 
values over all fully-resolved trees with n taxa. Fig. 3 plots the range of Zt for the values of /3 
used in Fig. 2, and for n — 10, 20, 30, 40, 50 taxa trees, on a log-log axis. The trees minimising Zt 
were always caterpillar trees and the trees maximising Zt were usually, but not always, balanced 
trees. The figure indicates that when f3 is outside the range [0.03,3] there is little variation in 
Zt between different trees. With 50 taxa, the normalising constants differ by a maximum of 7.5 
log- units. 



12 



DAVID BRYANT AND MIKE STEEL 



Suppose that we are comparing the log-likelihood of two trees Ti and T2 with respect to a third 
tree T. If dRF{T,Ti) ^ dRF{T,T2) then 

I log(e-'5<'(^'^i)) - log(e-^'^(^'^^))| > 2/3 

so ignoring the normalising constant will only change the order of likeUhood values if | log Zt^ — 
log^T2l ^ 2/3. Plotting the curve for 2/3 on Fig. 3 wc see that | log Zt^ — log .Era | > 2/3 for some 
pairs of 50-taxa trees only when /3 lies in the interval [1.25,1.86]. The corresponding interval 
will be even smaller for trees with fewer taxa: for 20 taxa trees there is no value of /3 for which 
ignoring Zt scores leads to a switch in the order of likelihood values for two trees. 

In summary, when f3 is approximately 1.5, and the number of taxa is greater than around 20, 
it is potentially important to correctly compute normalisation constants. Outside that range, 
the influence of Zt on likelihood rankings can be safely ignored. Wc note, however, that here 
we are only interested in relative ordering of supertrees with respect to likelihood: a Bayesian 
Monte-Carlo approach may well need accurate Zt values for all /3. 
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Figure 1. Two fully resolved phylogenetic trees on six leaves, with Robinson- 
Foulds distance two. 
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Figure 2. The average Zt values for different values of /3, plotted on a log- 
log axis. The approximations for small and large (3 are plotted. All values 
were computed by drawing 1000 fifty taxa trees from a uniform distribution and 
computing normalising constants exactly using the algorithms described here. 
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Figure 3. The range of the Zt values computed for different /3 and pfotted 
on a log-log axis. The Zt values were computed by drawing 1000 trees from 
a uniform distribution with n — 10,20,30,40,50 taxa (five curves). The range 
is the difference between the maximum Zt and minimum Zt values, for each 
choice of (3 and n. The dotted line indicates the 2(3 value: when the range is less 
than 2/3 ignoring the normalising constant has no effect on the relative order of 
likelihood values. 
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